System and method for admittance model identification for inverter-based resources

ABSTRACT

A system and method for obtaining frequency-domain admittance/impedance using a few sets of time-series data. Utilizing an Eigensystem Realization Algorithm (ERA) or a dynamic mode decomposition (DMD), the input/output frequency-domain model reflecting voltage/current relationship is identified, and admittance/impedance identification is demonstrated for a grid-connected inverter-based resource (IBR) system. The proposed approach provides a powerful tool to replace the state-of-the-art harmonic injection methodology.

BACKGROUND OF THE INVENTION

Unprecedented dynamic phenomena appear in power grids due to integrationof more and more inverter-based resources (IBR) (e.g., wind and solar) Amajor challenge is that inverter models are proprietary information andusually only real code models are provided to grid operators. Thus,measurement-based characterization of IBRs is a popular approach to findthe frequency domain measurements of an IBR's admittance or impedance.The predominant methods rely on injecting perturbation and extractingfrequency-domain information via Fast Fourier Transformer (FFT).

High penetration of IBRs introduce many new dynamic phenomena into powergrids. First, severe wind farm subsynchronous oscillations (SSO) due towind interactions with RLC circuits occurred in Texas in 2009. Similarphenomena were observed in North China in the years that followed. In2017, three more events were observed. Oscillations were also observedin real-world wind farms with weak grid interconnection. For voltagesource converter (VSC)-based HVDC (high-voltage direct current) withweak grid interconnection, stability issues have also been identified.

While electromagnetic (EMT)-type simulation in environments such asPSCAD and MATLAB/SimPowerSystems, is the major tool for dynamic study,EMT simulation mainly serves as a validation tool, instead of ananalysis tool. To facilitate small-signal analysis, linear models aredesired.

Details of inverters are usually proprietary information and are notcompletely released to grid operators. Measured impedances, throughknown frequency scanning or harmonic injection methods, thus serve asimportant characteristics for those devices. Impedance-based stabilityanalysis has been popularly used in the power electronics field as wellas type-3 wind SSO screening studies.

Currently, the predominant technology for ISR characterization isfrequency-domain measurement. Two types of methods are employed to findimpedance/admittance measurements. The first type, the harmonicinjection method or the frequency scanning method relies on setting up ameasurement testbed for an IBR only. A device is first connected to a 60Hz voltage source to operate at certain operating condition. To havevoltage perturbation, an additional voltage source with small magnitudeat a certain frequency is connected in series with the original voltagesource. The current component at that frequency will be measured usingFast Fourier Transformation (FFT) and the current/voltage phasorrelationship at this frequency is then obtained. An additional step,vector fitting must be carried out to obtain an input/output modeldescribed in either state-space or transfer function.

The second type does not require a measurement testbed. Rather, anadditional hardware device, impedance measurement unit (IMU), isinstalled in a system in operation. It injects series voltageperturbation, collects both voltage and current measurements'frequency-domain information through FFT and computes dq-frame impedancemeasurements.

Both methods require repeated sinusoidal injections and FFT analysis. Tospeed up the process, different types of injection signals have beenproposed and implemented in IMUs, e.g., chirp signal. Nevertheless, thecurrent IMUs all rely on injections and FFT to lead to impedance atfrequency points.

Accordingly, what is needed in the art is an improved system and methodfor generating an impedance or admittance model of an inverter-basedresource (IBR) that does not rely on numerous experiments requiringharmonic injection and frequency scan methods.

SUMMARY OF THE INVENTION

In various embodiments, the present invention provides a system andmethod for time-domain measurement-based admittance identification sothat event data can be utilized for admittance identification, inreal-time.

In one embodiment, the present invention provides a method foradmittance identification of a grid-connected inverter-based resource(IBR), which includes, capturing transient data in response to at leasttwo independent time-domain events on a bus between a grid-connected IBRand a power grid and identifying a dq-frame admittance model of the IBRfrom the captured transient data.

In a particular embodiment, the at least two independent time-domainevents are voltage perturbations on the bus and the method includescapturing and recording a transient current flowing to the IBR and atransient voltage time-series data at a terminal voltage of the IBR.

The method may further include, applying an Eigensystem RealizationAlgorithm (ERA) to the sampled transient data to generate an s-domainexpression of the transient data. Alternatively, the method may include,applying dynamic mode decomposition (DMD) to the sampled transient datato generate an s-domain expression of the transient data.

In an additional embodiment, the invention provides a system formeasuring an admittance of a grid-connected inverter-based resource(IBR). The system includes a measurement unit coupled to a point ofinterconnection (POI) between a power grid and an inverter-basedresource (IBR). The measurement unit is configured for capturingtransient data in response to at least two independent time-domainevents on a bus between the IBR and the power grid and for identifying adq-frame admittance model of the IBR from the captured transient data.

The measurement unit may comprise hardware circuitry and softwarecomponents that are well known in the art for measuring voltages andcurrents at a point of interest.

The embodiments of the present invention replace the harmonic injectionand frequency scan methods currently known in the art that requirenumerous experiments. The present invention thereby provides an improvedsystem and method for performing admittance measurements forinverter-based resources (IBRs), such as wind and solar farms.

BRIEF DESCRIPTION OF THE DRAWINGS

For a fuller understanding of the invention, reference should be made tothe following detailed disclosure, taken in connection with theaccompanying drawings, in which:

FIG. 1 is a block diagram of an analytical dq-frame model, in accordancewith an embodiment of the present invention.

FIG. 2 is a block diagram schematically depicting a measurement testbedfor a grid-connected voltage source converter (VSC) circuit, inaccordance with an embodiment of the present invention.

FIG. 3 illustrates Bode diagrams providing a graphical illustration ofthe admittance model frequency-domain responses based upon theanalytical model shown in FIG. 1.

FIG. 4 is a graphical illustration of the dynamic responses of i_(d) andi_(q) for two events. Event 1: v_(gd) step change of 0:001 pu at t=1 s.Event 2: v_(gq) step change of 0:001 pu at t=1 s, in accordance with anembodiment of the present invention.

FIG. 5A is a graphical illustration of actual eigenvalues of the 9thorder system and the estimation from current signals, in accordance withan embodiment of the present invention.

FIG. 5B is a zoomed-in view of the graphical illustration of FIG. 5A.

FIG. 6 is a graphical illustration of dynamic responses of the VSC gridintegration system subject to 0.5% chance in dc-link voltage order,illustrating that the system is unstable with growing 5 Hz oscillations,in accordance with an embodiment of the present invention.

FIG. 7 is a graphical comparison of eigenvalues obtained using anadmittance-based approach and an analytical model-based approach, inaccordance with an embodiment of the present invention.

FIG. 8A illustrates an electromagnetic transient (EMT) testbedcomprising a 400-kW grid-interconnected photovoltaic (PV) farm, inaccordance with an exemplary embodiment of the present invention.

FIG. 8B illustrates an electromagnetic transient (EMT) testbedcomprising a 100-MW grid-connected type-4 wind farm, in accordance withan exemplary embodiment of the present invention.

FIG. 9A is a graphical illustration of 7.5 Hz oscillation observed inthe EMT testbed of FIG. 8A on the AC side.

FIG. 9B is a graphical illustration of 7.5 Hz oscillation observed inthe EMT testbed of FIG. 8A on the DC side.

FIG. 10 is a graphical illustration of two sets of data generated byapplying small signal perturbations to the EMT testbed of FIG. 8A.

FIG. 11 illustrates Bode diagrams illustrating a comparison of 50^(th)order model and 10^(th) order model from dynamic mode decompositions(DMD), in accordance with an embodiment of the present invention.

FIG. 12A is a graphical illustration of reconstructed signals verses theoriginal measurements, wherein the dotted line represents the originalmeasurements, the thin solid line represents the reconstructed signals.

FIG. 12B is a zoomed-in version of FIG. 12A.

FIG. 13 is a graphical illustration comparing the Bode plots of theadmittance models obtained from ERA, DMD and the harmonic injectionmethod.

FIG. 14 is a graphical illustration of eigenvalues of the closed-loopsystem based on the 10^(th) order system from DMD, in accordance with anembodiment of the present invention.

FIG. 15A is a graphical illustration of event data resulting from aresistive load disconnection, in accordance with an embodiment of thepresent invention.

FIG. 15B is a graphical illustration of event data resulting from acapacitive load disconnection, in accordance with an embodiment of thepresent invention.

FIG. 16 is a graphical illustration of the dq-frame admittancemeasurements of the wind farm, the line and the true admittance of theline, in accordance with an embodiment of the present invention.

DETAILED DESCRIPTION OF THE INVENTION

Modern renewable energy resources generating devices (e.g. solar andwind), are interfaced to the grid through inverter technology, such asinverter-based resources (IBR). Inverters and IBRs consist of powerelectronics that drive the electrical control and performance of theseresources. Measured impedances of the inverter-based resources serve asimportant characteristics for impedance-based stability analysis. Theresults of the impedance-based stability analysis can be used to controlthe operation of the IBR and improve the grid stability.

The state-of-the-art technology for obtaining an impedance model of adevice requires numerous repeated harmonic injections to generate a Bodeplot. To overcome the deficiencies of the state-of-the art technology,in various embodiments, the present invention provides a system andmethod for obtaining a dq-frame admittance from time-domain signals.Dq-frame admittance has been shown to be able to accurately predictstability. In addition, frequency coupling phenomenon observed forconverters in the static frame is not a concern in the dq-frame.

Compared to a known frequency scanning method, which relies on ameasurement testbed and which requires repeated sinusoidal injections,leading to measurements at frequency points, a simpler experiment designis desired. Further, an input/output model, such as an s-domain model isdesired. With the s-domain admittance model, eigenvalue analysis thatpredicts system stability is possible.

Compared to another known method that requires an additional hardwaredevice to inject perturbations to a system to obtain measurements, thereare plenty of event data generated during IBR operation in a power grid,which can be utilized by the present invention to obtain an admittancemodel. A method of admittance identification without any perturbationinjection and capable of utilizing transient response data captured forevents is appealing.

In the present invention, ERA (eigensystem realization algorithm) andDMD (dynamic mode decomposition) are used to convert time-domain data toa corresponding s-domain expression. This capability further opens thedoor to an application with a high practical value: using event data toextract an MR's admittance in real-time.

In various embodiments, the present invention provides a novel procedurethat leads to dq-frame admittance identification using time-domain eventdata, e.g., step response data or load tripping data. In this procedure,two sets of data from two events lead to a dq-frame admittance. Themethod relies on a critical technology that converts time-domain signalsinto s-domain expressions. This technology employs ERA or DMD toaccurately estimate eigenvalues and residuals to construct the s-domainexpressions. In addition, since IBRs employ power electronic converters,and harmonics are common in signals, an efficient denoising technique tocancel harmonics has been proposed, tested and shown efficacy.

The present invention additionally provides a real-time dq-frameadmittance identification method using data from multiple events.Compared to the existing impedance measurement (IMU) technology, theproposed method does not require any additional hardware device toinject perturbations.

ERA assumes that the dynamic response is due to and impulse input.Consider a Linear-Time Invariant (LTI) system in the discrete domain asthe following:x _(k+1) =Ax _(k) +Bu _(k) ,y _(k) =Cx _(k) +Du _(x)  (1)

where y∈

^(K×1) is defined as the output column vector of the system with Koutput channels, A∈

^(n×n), B∈

^(n×1), C∈

^(K×n), and D∈

^(K×1) are system matrices.

It is noted that ERA does not lead to an input/output model due to theassumption of impulse input. Rather, ERA leads only to the s-domainexpression of the measurement data. Input to ERA is the measurement datasampled at equal intervals, notated as y₁, y₂, . . . y_(N). The outputfrom ERA is the state-space model defined by A, B, C, D matrices. Therank of the system n is an assumption fed to the algorithm. Eigenvaluesof the continuous system co can be found as ω_(j)=ln(λ_(j))/Δt, whereλ_(j) is an eigenvalue of A and Δt is the time interval.

The s-domain expression of y can be found as:

$\begin{matrix}{{y(s)} = {{{C( {{sI} - A} )}^{- 1}B} = {{\underset{\underset{\Phi}{︸}}{CV}( {{sI} - \Omega} )}^{- 1}\underset{\underset{b}{︸}}{V^{- 1}B}}}} & (2)\end{matrix}$

where V is the right eigenvector matrix of A and Ω is the diagonalmatrix of continuous system eigenvalues.

For k^(th) measurement, its s-domain expression can be found as:

$\begin{matrix}{{y_{k}(s)} = {\sum\limits_{j = 1}^{n}\frac{\phi_{kj}b_{j}}{s - \omega_{j}}}} & (3)\end{matrix}$

Assuming x₀=0, the system response due to an impulse input (u₀=1,u_(k)=0, k>0) can be found as follows:x ₀=0,x ₁ =B,x ₂ =AB, . . . ,x _(k) =A ^(k-1) By ₀ =D,y ₁ =CB,y ₂ =CAB, . . . ,y _(k) =CA ^(k-1) B

A critical step of ERA is to construct two shifted Hankel matrices.

$\begin{matrix}{H_{1} = \begin{bmatrix}y_{1} & y_{2} & \ldots & y_{L} \\y_{2} & y_{3} & \ldots & y_{L + 1} \\\vdots & \vdots & \ddots & \vdots \\y_{N - L + 1} & y_{N - L + 2} & \ldots & y_{N}\end{bmatrix}_{{K{({N - L + 1})}} \times L}} & (4) \\{H_{2} = \begin{bmatrix}y_{2} & y_{3} & \ldots & y_{L + 1} \\y_{3} & y_{4} & \ldots & y_{L + 2} \\\vdots & \vdots & \ddots & \vdots \\y_{N - L + 2} & y_{N - L + 3} & \ldots & y_{N}\end{bmatrix}_{{K{({N - L + 1})}} \times L}} & (5)\end{matrix}$

It can be seen that the Hankel matrices can be decomposed as follows:

$\begin{matrix}{H_{1} = {\begin{bmatrix}{CB} & {CAB} & \ldots & {CA^{L - 1}B} \\{CAB} & {CA^{2}B} & \ldots & {CA^{L}B} \\\vdots & \vdots & \vdots & \vdots \\{CA^{L - N}B} & {CA^{L - N + 1}B} & \; & {CA^{N - 1}B}\end{bmatrix} = {\underset{\underset{\mathcal{O}}{︸}}{\begin{bmatrix}C \\{CA} \\\vdots \\{CA^{N - L}}\end{bmatrix}}\underset{\underset{\mathcal{C}}{︸}}{\begin{bmatrix}B & {AB} & \ldots & {A^{L - 1}B}\end{bmatrix}}}}} & (6) \\{\mspace{79mu}{ \Rightarrow H_{1}  = {\mathcal{O}\mathcal{C}}}} & (7) \\{\mspace{79mu}{ \Rightarrow H_{2}  = {\mathcal{O}\; A\;\mathcal{C}}}} & (8)\end{matrix}$

where

is the observability matrix and

is the controllability matrix. ERA employs singular value decomposition(SVD) and rank reduction to find two matrices to realize

and

. Further, A matrix can be found by exploring the relationship betweenH₁ and H₂. Details of ERA implementation are well known in the art andcan be found in the literature, and as such, have been omitted from thisdisclosure.

Though ERA leads to a state matrix, the state variables are unknown. Onthe other hand, DMD results in a state matrix A(x_(k+1)=Ax_(k)) with thestate variables x clearly defined as the measurements, where x_(k)∈

^(n) and A∈

^(n×n). The time interval between two consecutive snapshots is Δ_(t).

Applying eigen-decomposition to matrix A(A=ΦΛΦ⁻¹) leads to:

$\begin{matrix}{x_{k + 1} = {{\Phi\Lambda\Phi^{- 1}x_{k}} = {{\Phi\Lambda^{2}\Phi^{- 1}x_{k - 1}} = {\ldots = {{\Phi\Lambda}^{k}\underset{\underset{b}{︸}}{\Phi^{- 1}x_{1}}}}}}} & (9)\end{matrix}$

where b∈

^(n×1), Φ∈

^(n×n) is the right eigen vector matrix of A and A is a diagonal matrixwith elements as λ_(i), i=1, . . . , n.

Equation (9) can also be written as follows:

$\begin{matrix}{x_{k + 1} = {\sum\limits_{j = 1}^{n}{\phi_{j}\lambda_{j}^{k}b_{j}}}} & (10)\end{matrix}$

The time-domain expression of x(t) can be constructed using theeigenvalues and eigenvectors of A.

$\begin{matrix}{{x(t)} = {{\sum\limits_{j = 1}^{n}{\phi_{j}e^{\omega_{j}t}b_{j}}} = {\Phi e^{\Omega t}b}}} & (11)\end{matrix}$

where Ω is a diagonal matrix that contains the continuous-eigenvalues,ω_(j). The relationship between the discrete and continuous eigenvalueis ω_(j)=ln(λ_(j))/Δt.

The inputs to the DMD algorithm are measurement data sampled at equalintervals, denoted as x₁, x₂, . . . , x_(N). The outputs are: theeigenvector matrix 1, the eigen value matrix Ω, and the initial stateprojected to the eigenvector basis b. Rank of the system r should alsobe assumed.

A critical step in DMD is to form a data matrix X. This matrix isconstructed to contain state for N snapshots with equal time-intervals.

$\begin{matrix}{X = \begin{bmatrix}| & | & \; & | \\x_{1} & x_{2} & \ldots & x_{N} \\| & | & \; & |\end{bmatrix}} & (12)\end{matrix}$

Following the construction of the data matrix, two shifted matrices areformed:

$\begin{matrix}{{X_{1}^{N - 1} = \begin{bmatrix}| & | & \; & | \\x_{1} & x_{2} & \ldots & x_{N - 1} \\| & | & \; & |\end{bmatrix}},} & ( {13a} ) \\{{X_{2}^{N} = \begin{bmatrix}| & | & \; & | \\x_{3} & x_{3} & \ldots & x_{N} \\| & | & \; & |\end{bmatrix}},} & ( {13b} )\end{matrix}$where X₁ ^(N−1)∈

^(n×(N−1)), X₂ ^(N)∈

^(n×(N−1)). The subscript and superscript refer to the first and lastmeasurement snapshots in the set, respectively. It can be seen that:X ₂ ^(N) =AX ₁ ^(N−1)

Matrix A can then be found by exploring the above relationship. Furtherdetails are omitted.

With DMD's outputs Φ, b and the eigenvalue vector co for continuoussystem, it is easy to find the s-domain expression of kth measurementsignal as:

${x_{k}(s)} = {\sum\limits_{j = 1}^{r}\frac{\phi_{kj}b_{j}}{s - \omega_{j}}}$

To illustrate proof-of-concept of the present invention, data isgenerated from an analytical model build in dq-frame. With theidentified admittance, stability analysis is also illustrated for toverify the design.

The dq-frame admittance of a grid-connected voltage source controller(VSC) 200 shown in FIG. 2 will be described. As illustrated, the VSC 200is coupled to a photovoltaic or wind resource 205 and to a power grid215 at a point of common coupling (PCC) 210, also known as a point ofinterconnection (POI)A VSC 200 with a grid-following control structureis assumed in this exemplary embodiment. The dq-frame model shown inFIG. 1 was designed and implemented. A small-signal model can beextracted from the analytical model at an operating condition throughnumerical perturbation. With the dq-frame grid voltage designated as aninput and the current flowing out of the inverter as the output, theinput/output model extracted by numerical perturbation is indeed −Y,where Y notates the admittance model. Thus, the admittance model forsuch a system is known, as shown in FIG. 3, where the admittance modelfrequency-domain responses are illustrated. This model can serve as thebenchmark for the identified admittance.

The foremost challenge is how to create data that is suitable foradmittance identification. In the harmonic injection method, sinusoidalperturbation is used. On the other hand, for time-domain data, Laplacetransform of the step response of a system is associated with theproduct of the transfer function of the system and 1/s. Thus, stepchanges may be used as perturbations.

A step change with 0.001 pu size will be applied to v_(gd). Linecurrents are measured and notated as i_(d) ⁽¹⁾ and i_(q) ⁽¹⁾. Anotherstep change with 0.001 pu size is applied to v_(gq) and the linecurrents are notated as i_(d) ⁽²⁾ and i_(q) ⁽²⁾. The step response datais presented in FIG. 4.

Data in the time frame from 1 second to 1.98 is used for analysis.Sampling frequency is 2500 Hz. Since ERA assumes the initial statevariables are zeros, the data are processed to have the initialsteady-state values taken off. Proper scaling is applied to signals tohave similar degree of variation. For the four signals, the scales are1000, 1, 100, and 10. The two sets of scaled event data are fed into theERA.

The estimated system is assumed to have an order of 10. This order ischosen by considering the 9th order system and the step responseperturbation. Given the step change input, the measurement signal shouldhave an order of 10.

The system matrix A is first computed based on the measurement signals.The eigenvalues of the continuous system are estimated and the residualsfor each signal can be computed. With the eigenvalue and residualsfound, measurement data can be reconstructed, and the transfer functioncan be found for each signal. FIG. 4 presents the reconstructed signalswith their initial steady-state values added. The reconstructed signalshave a close to 100% match with the original data.

FIG. 5 presents the estimated eigenvalues versus the eigenvalues of theoriginal system in the complex plane. It can be seen that the identifiedeigenvalues align very well with the original eigenvalues. The estimatedeigenvalue at the original point is due to the step input whichintroduces an “s” in the denominator of the current transfer functions.

For each signal, its estimated Laplace domain expression can now befound: i_(d) ⁽¹⁾(s), i_(q) ⁽¹⁾(s), i_(d) ⁽²⁾(s) and i_(q) ⁽²⁾(s), wheresuperscript notates event number. The admittance model can be found bytaking into the effect of step response of v_(d) or v_(q).

$\begin{matrix}{Y = {\frac{- s}{p}\begin{bmatrix}{i_{d}^{(1)}(s)} & {i_{d}^{(2)}(s)} \\{i_{q}^{(1)}(s)} & {i_{q}^{(2)}(s)}\end{bmatrix}}} & (14)\end{matrix}$

where p is the size of perturbation. For this study, p=0:001.

FIG. 3 presents the frequency responses of the estimated admittanceversus the original admittance. It can be seen that the estimatedadmittance provides an excellent match in the frequency range of 1 Hz to100 Hz. Considering that grid dynamic studies are concerned aboutdynamics in this range, the proposed approach, using step responses of 1second data for admittance identification, has been demonstrated as apowerful tool for device characterization.

The analytical model shown in FIG. 1 can be used to demonstratelow-frequency oscillations when the VSC is integrated into a weak gridat X_(g)=0.8 pu. FIG. 6 presents the system response subject to a smallstep change in dc-link voltage order. It can be seen that the system issubject to 5 Hz oscillations.

The system is viewed at PCC (point of common coupling) bus with twoshunt admittances: Y_(VSC) and Y_(line). Applying circuit analysis, itcan be found that the relationship between the injected small current atthe PCC bus and the PCC bus deviation is as follows:

$\begin{matrix}{\begin{bmatrix}{\Delta i_{d,{inj}}} \\{\Delta i_{q,{inj}}}\end{bmatrix} = {( {Y_{VSC} + Y_{line}} )\begin{bmatrix}{\Delta V_{d,{PCC}}} \\{\Delta V_{q,{PCC}}}\end{bmatrix}}} & (15)\end{matrix}$

If (15) is viewed as an input/output system, then the PCC voltage is theoutput (notated as y) and the injected current is the input (notated asu). Hence, the transfer function matrix notates as G(s) from u to y is:

$\begin{matrix}{{G(s)} = {\frac{y}{u} = \underset{\underset{Y}{︸}}{( {Y_{VSC} + Y_{line}} )^{- 1}}}} & (16) \\{{{where}\mspace{14mu} Y_{line}} = \begin{bmatrix}{R_{g} + {s\; L_{g}}} & {{- \omega_{0}}L_{g}} \\{\omega_{0}L_{g}} & {\;{R_{g} + {s\; L_{g}}}}\end{bmatrix}^{- 1}} & (17)\end{matrix}$

and ω₀ is the nominal frequency 377 rad/s.

Poles of G(s) are the eigenvalues of the system. In turn, roots ofdet(Y) or the zeroes of the s-domain admittance matrix Y are theeigenvalues of the system. With the VSC admittance identified throughtime-domain signals and the line admittance known, eigenvalues of thesystem can be found.

FIG. 7 presents the eigenvalues of the admittance-based approach and theoriginal set. As shown, the figure illustrates an exact match.

While the data generated from the analytical model contain little noise,data generated at the real world contain noise and harmonics. Inaddition, the real-world physical system is usually of high order. Thus,a more sophisticated electromagnetic transient (EMT) testbed (EMTTestbed 1) is presented below for demonstration. A notable procedure ondenoising is described for the ERA and DMD identification tools.

In general, the present invention provides a measurement unit that canbe connected at the point of interconnection (POI) between a power gridand a solar farm or wind farm. The measurement unit captures and recordstransient data of voltage at the POI and the current flowing into theinverter-based resource (IBR). With voltage and current data from atleast two time-domain events, an admittance model of the IBR (solar farmor wind farm) can be found.

From the admittance model, one can quickly project if the IBRinterconnected system is stable or not at a certain power grid strengthcondition. The assumption is that the admittance model of the point ofinterconnection is know if a specific grid strength is assumed. One canuse the admittance-based stability check for this speculation. Equations(15)-(17) illustrate how to obtain aggregated admittance of the entiresystem after the IBR admittance has been measured with the measurementunit. Knowing the aggregated admittance, eigenvalues of the system canbe found, and a stability check can be made, as demonstrated in FIG. 5and FIG. 14.

For example, FIG. 14 indicates that for the 400-kW solar farm shown inFIG. 8A, once its admittance is measured, one can project that it canwork in a grid with Thevenin equivalent impedance less than 1.02. If thegrid is weaker, or the impedance is greater than 1.02, then the systemloses stability and 7.7 Hz oscillations will appear.

The method of the present invention for deriving the admittance model ofthe IBR can be applied to any subsystem, as long as the time-domainevent that trigger the transients occur outside of the subsystem.

EMT Testbed 1, shown in FIG. 8A, is build based on the demo 400-kW PVsystem in MATLAB/SimPowerSystems. This testbed has full dynamicsmodeled, including not only the grid-connected VSC average model and itscontrol, transmission line electromagnetic dynamics, but also dc/dcboost converter average model, dc/dc boost converter control, maximumpower point tracking (MPPT) and dc-side circuit dynamics. The PV voltage(V_(PV)) at the input side of the dc/dc boost converter is 260V. Outputof the dc/dc boost converter is at 500V. The 500V dc voltage isconverted to a three-phase 260 V ac power voltage by a VSC. A shuntcapacitor (C₁) is connected on the point of common coupling (PCC) bus tocompensate reactive power. There are two Y ground transformers (T₁, T₂)on transmission line and the system is connected to a 120 kV grid bus.Note that R_(g) and L_(g) are the resistance and inductance on gridtransmission line. The system parameters are listed in Table I.

TABLE 1 PV System Parameters Description Parameters Value Power baseS_(base) 400 kW Power level P 0.935 pu System frequency f_(base) 60 HzConverter filter R₁ 0.15/50 pu X₁ 0.15 pu Shunt capacitor C₁ 0.25 puDC-link capacitor C_(dc) 0.054 F Transmission line R_(g) 0.1 X_(g) Innerloop K_(pi), K_(ii) 0.3, 5   DC-link control loop K_(pp), K_(ip)  1, 100V_(ac) control loop K_(pv), K_(iv)  1, 100 PLL K_(p,PLL), K_(i,PLL)  60,1400

7.5 Hz low-frequency oscillations can be observed in FIG. 9 when thetransmission line impedance X_(g) is 1.0 pu and the PCC voltage controlorder steps from 1.0 pu to 0.99 pu at t=4 s. Both the ac side and the dcside measurements are presented in FIG. 9.

The admittance measurement testbed is set up to have the PV farmconnected to a small impedance (X_(g)=0.1 and R_(g)=0.01 pu) to avoltage source. The voltage source's v_(gd) and v_(gq) are configured sothat the real power, reactive power from the PCC bus to the grid remainsthe same between the EMT testbed and the measurement testbed.

Two sets of data are generated by applying small signal perturbations(0.02 pu) on v_(g,d) ^(g) and v_(g,q) ^(g), respectively. FIG. 10presents the data. The data are resampled to having a sampling frequencyof 4 kHz. It can be observed that noises and harmonics exist even atsteady-state.

ERA and DMD tools will be applied to handle the two sets of the data(each set has two time-domain signals) to identify the s-domainadmittance and reconstruct the signals.

For DMD, the measurement data from 1 second to 1.5 seconds, with theirinitial values taken off, will be used. As a total, there are fourmeasurements. The data matrix X will have a dimension of 4×2000. On theother hand, the order of the system is much higher than 4 and DMD willnot be able to lead to a system matrix of higher order.

In this exemplary embodiment, the stacking number s_(t) is chosen to behalf of the total snapshot number (N=2000 and s_(t)=1000). This numberis also used in the ERA exemplary embodiment. The dimension of the datamatrix becomes 4000×1000. Initially, 50-order rank is assumed. DMDoutputs a vector b which indicates each eigenvalue's role in the dynamicsystem. A small absolute value indicates negligible role. Hence a filtercan be set up to filter out those eigenvalues whose |b_(j)| is less than1% of the maximum 1% max(|b|). The technique is comparable to thedenoising technique relying on FFT coefficients. With this filteringprocedure, a 10th order model is obtained.

FIG. 11 presents the Bode plots of the 50th order model and the 10thorder model. It can be seen that the reduced order model leads to muchsmoother Bode plots.

The measurement signals can be reconstructed using the reduced-ordersystem and various methods known in the art. FIG. 12 presents theoriginal signals in dash-dotted lines and the reconstructed signals inblack thick solid lines. To provide a better view, a zoom-in diagram isalso provided. It can be observed that excellent match has beenachieved.

In another embodiment, ERA is also applied to the data set to carry outidentification. Similarly, a denoising procedure is also applied.

For the EMT testbed with instantaneous voltage and currents as statevariables, it is not possible to directly obtain a linearized model vianumerical perturbation since state variable are periodic. To validatethe admittance identified by DMD and ERA, the harmonic injection methodis applied to measure the admittance.

FIG. 13 presents the comparison of the two s-domain models with theadmittance measurement from the harmonic injection method. They showexcellent match.

It should be noted that the DMD and ERA identified models are s-domainmodels and can be directly used for eigenvalue computing relying ons-domain admittance. The harmonic injection leads to measurements. Toobtain s-domain model from the harmonic injection method, a vectorfitting procedure needs to be carried out.

Finally, stability analysis for weak grid operation is carried out andthe closed-loop system eigenvalues are shown in FIG. 14. It can be seenthat when the grid strength reduces, an oscillation mode at 7 Hz movesto the right half plane and X_(g)=1.02 is identified as a marginalcondition. This result corroborates with the EMT simulation results inFIG. 9.

In an exemplary embodiment, admittance measurement of a type-4 wind farmusing two sets of event data is presented. A 100-MW wind farm gridintegration system is shown in FIG. 8B. The EMT testbed is developedbased on the type-4 wind farm demo system in MATLAB/SimPowerSystems. Thetype-4 wind turbine model includes a synchronous generator, themechanical system and control, a three-phase rectifier, a dc-dc boosterand control, and a dc/ac converter and control. Average models areadopted for the dc/dc booster and the dc/ac converter. The type-4 windfarm consists of 50 2-MW wind turbines and the terminal ac voltage levelis 575 V. This farm is connected to a grid through a 575 V/25 kV step-uptransformer. The exporting power of the wind farm is 100 MW with thewind speed assumes to be 15 m/s. The measurement data collected are allscaled to pu values.

The power base is 111 MW and the wind farm is exporting level is 0.9 puor 100 MW.

Its PCC bus voltage is controlled at 1 pu or nominal level. Thedistribution line has a total reactance X_(g)=0.4 pu and a totalresistance of R_(g)=0.04 pu. At a bus notated for v, a resistive loadconsuming 0.1 pu real power or a capacitor injecting 0.1 pu reactivepower may be connected. Two events are simulated.

For the first event, the bus is supplying the resistive load at thebeginning. At t=5 seconds, the resistive load is disconnected. Thethree-phase bus voltage v, the current to the wind farm i₁, and thecurrent to the grid i₂ are all measured and converted to the dq-frame.

For the second event, the bus is connected for the capacitor at thebeginning. At t=5 seconds, the capacitor is disconnected. Using theaforementioned procedure, the event data in the dq-frame are obtained.The two sets of the event data are presented in FIG. 15. The simulationdata has a sampling rate of 50 μs.

The two sets of the data will be notated using superscript (i) to notatewith event. For each set of the data, the s-domain expression of the sixsignals will be learned using ERA. 1 second data from 5 seconds to 6seconds are used for learning. The data are resampled to have a samplingfrequency of 2000 Hz. The system order is assumed to be 25.

With s-domain expressions of all measurements known, the dq-frameadmittance is to be sought. When the wind farm is operating at the samecondition, its admittance will be kept the same. Hence it is true thatthe current and voltage for different events are related with the sameadmittance:

$\begin{matrix}{\underset{\underset{I{(s)}}{︸}}{\begin{bmatrix}{i_{d1}^{(1)}(s)} & {i_{d1}^{(2)}(s)} \\{i_{q}^{(1)}(s)} & {i_{q}^{(2)}(s)}\end{bmatrix}} = {{Y_{wind}(s)}\underset{\underset{V{(s)}}{︸}}{\begin{bmatrix}{v_{d1}^{(1)}(s)} & {v_{d1}^{(2)}(s)} \\{v_{q}^{(1)}(s)} & {v_{q}^{(2)}(s)}\end{bmatrix}}}} & (18)\end{matrix}$

Thus, the wind farm admittance at any frequency ω can be found using thefollowing equation.Y _(wind)(jω)=I(jω)V(jω)⁻¹  (19)

If there are m events and m>2, we may have the current and voltagemeasurement matrices as fat matrices, with a column dimension of m:I(jω), V(jω)∈

^(2×m). Thus, pseudo inverse or Moore-Penrose inverse of the voltagemeasurement matrix may be used.

Eq. (19) can be changed as follows:Y _(wind)(jω)=I(jω)V(jω)^(†)  (20)

where † refers to Moore-Penrose inverse.

From the current measurements of the transmission line, one may alsofind the admittance measurement of the transmission line. For the EMTtestbed, the transmission line is represented by a series RL componentand the parameters are given as R_(g)=0.04 pu and X_(g)=0.4 pu. Hence,its true dq-frame admittance can be found using (17).

One can compare the line admittance from the measurement data with thetrue line admittance. If they agree with each other, that means themethod leads to reasonable admittances for both line and the wind farm.

FIG. 16 presents the dq-frame admittance measurements of the wind farm,the line and the true admittance of the line. It can be seen that forthe line, the admittance measurement and the true admittance match verywell for frequency below 100 Hz. Thus, it is also concluded that thewind farm admittance measurement below 100 Hz can be used for analysis.

The proposed method leads to real-time admittance measurement usingmultiple event data. Compared to the current impedance measurement unittechnology, this method does not require any additional hardware toinject perturbations into a system.

The innovation of the present invention is admittance identificationusing time-domain data via a critical step: obtaining s-domainexpressions of measurements. Compared to the state-of-the art approach,i.e., harmonic injection method, the proposed approach of the presentinvention is much more efficient. The system and method of the presentinvention requires only two sets of data to lead to a dq-frameadmittance. The proof-of-concept has been demonstrated using ananalytical model. Two IBR grid-integration EMT testbeds were utilized todemonstrate the two admittance identification procedures. Specifically,an application of using event data for real-time admittance identifiedhas been demonstrated.

Hardware and Software Infrastructure Examples

The present invention may be embodied on various platforms. Thefollowing provides an antecedent basis for the information technologythat may be utilized to enable the invention.

Embodiments of the present invention may be implemented in hardware,firmware, software, or any combination thereof. Embodiments of thepresent invention may also be implemented as instructions stored on amachine-readable medium, which may be read and executed by one or moreprocessors. A machine-readable medium may include any mechanism forstoring or transmitting information in a form readable by a machine(e.g., a computing device). For example, a machine-readable medium mayinclude read only memory (ROM); random access memory (RAM); magneticdisk storage media; optical storage media; flash memory devices;electrical, optical, acoustical or other forms of propagated signals(e.g., carrier waves, infrared signals, digital signals, etc.), andothers. Further, firmware, software, routines, instructions may bedescribed herein as performing certain actions. However, it should beappreciated that such descriptions are merely for convenience and thatsuch actions in fact result from computing devices, processors,controllers, or other devices executing the firmware, software,routines, instructions, etc.

The machine-readable medium may be, for example, but not limited to, anelectronic, magnetic, optical, electromagnetic, infrared, orsemiconductor system, apparatus, or device, or any suitable combinationof the foregoing. More specific examples (a non-exhaustive list) of thecomputer readable storage medium would include the following: anelectrical connection having one or more wires, a portable computerdiskette, a hard disk, a random access memory (RAM), a read-only memory(ROM), an erasable programmable read-only memory (EPROM or Flashmemory), an optical fiber, a portable compact disc read-only memory(CD-ROM), an optical storage device, a magnetic storage device, or anysuitable combination of the foregoing. In the context of this document,a computer readable storage medium may be any non-transitory, tangiblemedium that can contain, or store a program for use by or in connectionwith an instruction execution system, apparatus, or device.

A machine-readable signal medium may include a propagated data signalwith machine-readable program code embodied therein, for example, inbaseband or as part of a carrier wave. Such a propagated signal may takeany of a variety of forms, including, but not limited to,electro-magnetic, optical, or any suitable combination thereof. Amachine-readable signal medium may be any machine-readable medium thatis not a computer readable storage medium and that can communicate,propagate, or transport a program for use by or in connection with aninstruction execution system, apparatus, or device. However, asindicated above, due to circuit statutory subject matter restrictions,claims to this invention as a software product are those embodied in anon-transitory software medium such as a computer hard drive, flash-RAM,optical disk or the like.

Program code embodied on a machine-readable medium may be transmittedusing any appropriate medium, including but not limited to wireless,wire-line, optical fiber cable, radio frequency, etc., or any suitablecombination of the foregoing. Machine-readable program code for carryingout operations for aspects of the present invention may be written inany combination of one or more programming languages, including anobject oriented programming language such as Java, C#, C++, Visual Basicor the like and conventional procedural programming languages, such asthe “C” programming language or similar programming languages.

Aspects of the present invention are described below with reference toflowchart illustrations and/or block diagrams of methods, apparatus(systems) and computer program products according to embodiments of theinvention. It will be understood that each block of the flowchartillustrations and/or block diagrams, and combinations of blocks in theflowchart illustrations and/or block diagrams, can be implemented bymachine-readable program instructions.

The advantages set forth above, and those made apparent from theforegoing disclosure, are efficiently attained. Since certain changesmay be made in the above construction without departing from the scopeof the invention, it is intended that all matters contained in theforegoing disclosure or shown in the accompanying drawings shall beinterpreted as illustrative and not in a limiting sense.

What is claimed is:
 1. A method for admittance identification of agrid-connected inverter-based resource (IBR), the method comprising:capturing transient data in response to at least two independenttime-domain events on a bus between a grid-connected IBR and a powergrid; sampling the captured transient data to generate sampled transientdata; generating an s-domain expression of the transient data from thesampled transient data; and identifying a dq-frame admittance model ofthe IBR from the s-domain expression of the transient data.
 2. Themethod of claim 1, wherein the at least two independent time-domainevents comprises a first time-domain event and a second time-domainevent and wherein the method further comprises: capturing a first set oftransient data in response to the first time-domain event; and capturinga second set of transient data in response to the second time-domainevent.
 3. The method of claim 1, wherein the at least two independenttime-domain events are step change events.
 4. The method of claim 1,wherein the at least two independent time-domain events are voltageperturbations on the bus and wherein capturing the transient data inresponse to the at least two independent time-domain events on the buscomprises measuring a transient current flowing to the IBR.
 5. Themethod of claim 1, wherein the at least two independent time-domainevents are voltage perturbations on the bus and wherein capturing thetransient data in response to the at least two independent time-domainevents on the bus comprises measuring a transient voltage time-seriesdata at a terminal voltage of the IBR.
 6. The method of claim 1, whereinsampling the captured transient data to generate the sampled transientdata comprises sampling the captured transient data at equal intervalsand wherein generating the s-domain expression of the transient datafurther comprises applying an Eigensystem Realization Algorithm (ERA) tothe sampled transient data to generate the s-domain expression of thetransient data.
 7. The method of claim 1, wherein sampling the capturedtransient data to generate the sampled transient data comprises samplingthe captured transient data at equal intervals and wherein generatingthe s-domain expression of the transient data further comprises applyingdynamic mode decomposition (DMD) to the sampled transient data togenerate the s-domain expression of the transient data.
 8. The method ofclaim 7, further comprising denoising the s-domain expression of thetransient data.
 9. The method of claim 1, wherein the transient data iscaptured during real-time operation of the IBR.
 10. The method of claim1, wherein the at least two independent time-domain events are theresult of external perturbations outside of the IBR device.
 11. Themethod of claim 1, wherein the IBR is a wind resource or a photovoltaicsolar resource.
 12. A method for admittance identification of agrid-connected voltage source converter inverter-based resource (IBR),the method comprising: capturing a first set of transient data inresponse to a first time-domain event on a bus between a grid-connectedinverter-based resource (IBR) and) a power grid, wherein the first setof transient data comprises transient data for a first current flowingto the IBR in response to the first time-domain event on the bus and afirst transient data for a bus voltage; capturing a second set oftransient data in response to a second time-domain event on the busbetween the IBR and the power grid, wherein the second time-domain eventis independent of the first time-domain event and wherein the second setof transient data comprises transient data for a second current flowingto the IBR in response to the second time-domain event on the bus and asecond transient data for the bus voltage; and identifying a dq-frameadmittance model of the grid-connected IBR from the first set oftransient data and the second set of transient data.
 13. The method ofclaim 12, further comprising: sampling the captured transient data atequal intervals to generate sampled transient data; and applying anEigensystem Realization Algorithm (ERA) to the sampled transient data togenerate an s-domain expression of the transient data.
 14. The method ofclaim 12, further comprising: sampling the captured transient data atequal intervals to generate sampled transient data; and applying dynamicmode decomposition (DMD) to the sampled transient data to generate ans-domain expression of the transient data.
 15. The method of claim 13,further comprising denoising the s-domain expression of the transientdata.
 16. The method of claim 14, further comprising denoising thes-domain expression of the transient data.
 17. A system for measuring anadmittance of a grid-connected inverter-based resource (IBR), the systemcomprising: a measurement unit coupled to a point of interconnection(POI) between a power grid and an inverter-based resource (IBR); themeasurement unit for: capturing transient data in response to at leasttwo independent time-domain events on a bus between the IBR and thepower grid; sampling the captured transient data to generate sampledtransient data; generating an s-domain expression of the transient datafrom the sampled transient data; and identifying a dq-frame admittancemodel of the IBR from the s-domain expression of the transient data. 18.The system of claim 17, wherein the at least two independent time-domainevents are voltage perturbations on the bus and wherein capturing thetransient data in response to the at least two independent time-domainevents on the bus comprises measuring a transient current flowing to theIBR.
 19. The system of claim 17, wherein the at least two independenttime-domain events are voltage perturbations on the bus and whereincapturing the transient data in response to the at least two independenttime-domain events on the bus comprises measuring a transient voltagetime-series data at the POI.
 20. The system of claim 17, wherein the IBRis a wind resource or a photovoltaic solar resource.